///////////////////////////////////////////////////////////////////////////////
//
// nnMeteo.cc
//
// programma di esempio per l'uso e la verifica dei risultati di una
// rete neurale addestrata a fare previsioni metereologiche
//
// - converte un file dati contenente per un periodo di giorni consecutivi:
//   - pressione
//   - umidita'
//   - temperatura
//   - precipitazioni
//
// - definisce una rete neurale che usa i dati dei 4 gg precedenti per prevedere
//   un dato del giorno successivo (in questa versione, la temperatura)
//
// - i dati giornalieri vengono copiati in n-uple di Root in 5 file
//   - "dati.root" contiene i dati dal giorno no. 5 an no. N
//   - "dati-1.root" contiene i dati dal giorno n. 4 an no. N-1
//   - "dati-2.root" contiene i dati dal giorno n. 3 an no. N-2
//   - "dati-3.root" contiene i dati dal giorno n. 2 an no. N-3
//   - "dati-4.root" contiene i dati dal giorno n. 1 an no. N-4
//
// - viene costruita una n-upla estesa che contiene
//   - i dati del giorno per cui si vuole fare la previsione
//     (di questi dati verra' usato solo il dato che si vuole prevedere, sia
//      per addestrare la rete neurale sia per verificare l'addestramento)
//   - i dati dei giorni precedenti, usando la funzione Root "AddFriend"
//
// - viene fatto l'addestramento
// - vengono verificati i risultati su un campione indipendente di test
//
///////////////////////////////////////////////////////////////////////////////
//
// Istruzioni per l'uso:
//
// 1) eseguire root:
// shell> root -l <return>
//
// 2) compilare ed eseguire il programma da dentro Root
// root [0] .x nnMeteo.cc+
//
// (per eseguire in una volta sola la conversione dei dati (da ascii a
// root-tree), la visualizzazione delle variabili, e l'addestramento e
// utilizzo della rete.)
//
// 3) oppure, sempre da dentro Root:
//
// root [1] .L nnMeteo.cc+              # carica e compila il codece, necessario dopo ogni modifica
// root [2] converti()                  # converte dati.txt in in "dati.root", "dati-1.root", ...
// root [3] visualizza()                # visualizza la distribuzione dei dati in "dati.root" 
// root [4] reteneurale(algoritmo, neur1, neur2, cicli)
//   # addestra la rete usando le n-uple di Root in "dati*.root"
//   # - algoritmo	identifica l'algoritmo di addestramenti, da 1 a 6
//   # - neur1, neur2:  no. di neuroni nel primo e nel secondo strato interno (il 2o puo' essere zero)
//   # - cicli		numero di cicli di addestramento
// root [5] usarete()        (per utilizzare la rete gia' addestrata di cui avete salvato i pesi in pesi.txt, sul sample pippo.root)
//

//--- se non commentato, predici se piove o no
// #define PRED_PIOGGIA
//--- minima quantita' di pioggia per definire un giorno piovoso
#define PIOGGIA_MIN (1)

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <assert.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <numeric>

#include <string.h>

#include "TROOT.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TH1F.h"
#include "TMath.h"
#include "TFile.h"
#include "TTree.h"
#include "TEventList.h"
#include "TMultiLayerPerceptron.h"
#include "TMLPAnalyzer.h"

#endif

using namespace std;

//
// meteoData: classe per oggetto contenente i dati metereologici
//
// le macro ClassDef e ClassImp servono perche' Root possa
// usare le variabili e i metodi della classe per srivere/leggere
// su disco e per le formule [ es. ntupla->Draw() ]
//
// i dati si possono accedere sia con il vettore data[]
// sia con i metodi
// - pressione()
// - temperatura()
// - umidita()
// - precipitazione()
//
// le parentesi () sono necessarie ma a parte questo dettaglio questi
// metodi si usano come se fossero variabili, cioe' si possono usare
// per fare assegnamenti del tipo "md.pressione() = 1000;"
//
class meteoData {
 private:
 public:
  enum { meteoVarNum = 4 };
  Double_t data[meteoVarNum];

  //--- access for reading and writing
  Double_t& pressione() {return data[0];}
  Double_t& temperatura() {return data[1];}
  Double_t& umidita() {return data[2];}
  Double_t& precipitazione() {return data[3];}

  //--- access for reading when object is declared constant
  Double_t const& pressione() const {return data[0];}
  Double_t const& temperatura() const {return data[1];}
  Double_t const& umidita() const {return data[2];}
  Double_t const& precipitazione() const {return data[3];}

  meteoData() { memset( &data[0], sizeof(data), 0); }
  virtual ~meteoData() {}
  ClassDef(meteoData,1)  // per segnalare a Root il nome delle variabili
};

ClassImp(meteoData);

//
// salva i dati raccolti nel vettore dati in un file di Root
// 
//
//
void salvaNupla(TString name, vector<meteoData>& mdl, Int_t nSkip, Int_t nSkipMax )
{
  //--- il nome del file di Root dipende da quale giorno si parte
  TString output = name;
  //--- nome dell'oggetto nella n-upla
  TString mdName = "md";

  if (nSkip != 0) {
    output += "-";
    output += nSkip;
    mdName += nSkip;
  }
  output += ".root";
  
  // ora dobbiamo salvare tutto in un "tree" di root
  // apertura file di output
  TFile rfile(output,"RECREATE");

  TTree tree("T", "il nostro tree");
  meteoData* mdp = new meteoData();
  
  tree.Branch(mdName, "meteoData", &mdp, 64000, 1);
  
  //
  // produciamo un file contenente i dati corrispondenti ai giorni
  // da <nSkip> a <ndays - (nSKipMax-nSkip) -1> (inclusi)
  //
  Int_t ndays = mdl.size();
  Int_t iDay = 0;
  Int_t iValidDay = 0;
  
  vector<meteoData>::iterator iter;
  for( iter = mdl.begin(); iter != mdl.end(); iter++, iDay++ ) {
    //
    // salviamo nell'oggetto puntato da mdp
    // l'-i-esimo giorno di dati metereologici
    //
    *mdp = *iter;
    
    //
    // eliminiamo i dati non validi
    // questo significa che in alcuni casi il "giorno precedente"
    // corrispondera' in realta' a diversi giorni prima
    // non e' facile evitare questo senza significative
    // modifiche al programma
    //
    if (mdp->pressione() == -9999 ||
	mdp->temperatura() == -9999 ||
	mdp->umidita() == -9999 ||
	mdp->precipitazione() == -9999) {
      continue;
    }
    
    //
    // qui si possono regolarizzare i dati
    //

    //
    // se la quantita' di pioggia supera un valore minimo, si
    // definisce il giorno come "piovoso" con valore "1"
    // altrimenti il valore e' zero
    //
    mdp->precipitazione() = mdp->precipitazione() > PIOGGIA_MIN ? 1 : 0;
    
    //--- contatore per i soli giorni con dati validi
    if (iValidDay++ < nSkipMax -nSkip) continue;
    if (iDay >= ndays - nSkip) continue;

    //--- mettiamo i valori nella n-upla
    tree.Fill();
    // cout << "day " << iDay << " pressione " << mdp->pressione() << endl;
  }
  
  // operazioni conclusive: chiusura file di input, scrittura degli output
  tree.Write();
  rfile.Write();
  rfile.Close();
  delete mdp;
  
  cout << "Scritto il Root file " << output << endl;
}

//
// converte il file di testo name.txt nel file di Root name.root
// attenzione: se invocato piu' volte il file di Root risultante
// corrispondera' ai dati del solo ultimo file di dati
//
void converti(TString name = "dati", Int_t ndays=365)
{
  // definizione dei nomi dei file di input e output
  TString input = name + ".txt";

  // apertura file dati per leggere
  ifstream in;
  in.open(input);

  // lettura del contenuto del file di input
  vector<meteoData> mdl(ndays);
  meteoData md;
  
  Int_t nlines = 0;
  Double_t dummy;

  for(;;) {
    for (Int_t i=0; i<ndays; i++) {
      //
      // ogni linea corrisponde a 365 dati di una variabile
      // --> si aggiorna il campo corrispondente nell' i-esimo elemento
      // del vettore, la corrispondenza e':
      // linea 1 -> pressione
      // linea 2 -> temperatura
      // linea 3 -> umidita'
      // linea 4 -> precipitazione
      //
      if (nlines < meteoData::meteoVarNum) {
	in >> mdl[i].data[nlines];
      } else {
	// --- dati non utilizzati nelle linee successive a meteoVarNum
	in >> dummy;
      }
    }
    if (!in.good()) break;
    nlines++;
  }
  in.close();

  cout << "Letti dati corrispondenti a " << ndays << " giorni" << endl;

#if 0
  vector<meteoData>::iterator iter;
  for( iter = mdl.begin(); iter != mdl.end(); iter++) {
    cout
      << "p, t, h, r = "
      << iter->pressione() << " "
      << iter->temperatura() << " "
      << iter->umidita() << " "
      << iter->precipitazione() << " "
      << endl;
  }
#endif

  //
  // salva i dati in diversi file di Root contenenti gli stessi valori
  // ma sfasati da 0 a 4 giorni
  //
  salvaNupla(name, mdl, 0, 4);
  salvaNupla(name, mdl, 1, 4);
  salvaNupla(name, mdl, 2, 4);
  salvaNupla(name, mdl, 3, 4);
  salvaNupla(name, mdl, 4, 4);
  
  return;
}

void visualizza(TString name = "dati")
{
  TString rfname = name + ".root";

  // apertura del file e del suo tree
  TFile rfile(rfname);  
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)gDirectory->Get("T");

  TCanvas* c1 = new TCanvas("c1", "Variabili", 25,25, 700,700);
  c1->Divide(2,2);

  c1->cd(1);
  tree->Draw("pressione()");
  c1->cd(2);
  tree->Draw("temperatura()");
  c1->cd(3);
  tree->Draw("umidita()");
  c1->cd(4);
  tree->Draw("precipitazione()");
}

TString nomeAlgo(Int_t nAlgo)
{
  TString algo;
  if (nAlgo == 1) algo="Stochastic";
  if (nAlgo == 2) algo="Batch";
  if (nAlgo == 3) algo="SteepestDescent";
  if (nAlgo == 4) algo="RibierePolak";
  if (nAlgo == 5) algo="FletcherReeves";
  if (nAlgo == 6) algo="BFGS";

  return algo;
}

TString nomeBase(TString name,
		 Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain)
{
  // salva i pesi di questa rete per poterli usare dopo:
  TString base = name + "_";
  base += nomeAlgo(nAlgo);
  base += "_";
  base += neuro1;
  base += "_";
  base += neuro2;
  base += "_";
  base += ntrain;

  return base;
}

TString nomeFilepesi(TString name,
		     Int_t nAlgo, Int_t neuro1, Int_t neuro2,Int_t ntrain)
{
  return
    nomeBase(name, nAlgo, neuro1, neuro2, ntrain)
    + "_pesi.txt";
}

TString nomeNnConfig(TString name,
		     Int_t nAlgo, Int_t neuro1, Int_t neuro2,Int_t ntrain)
{
  return
    nomeBase(name, nAlgo, neuro1, neuro2, ntrain)
    + "_config.txt";
}

//
// costruisci la stringa che definisce la rete neurale
//
TString nnStruct(Int_t neuro1, Int_t neuro2)
{
  //
  // le istruzioni sono su http://root.cern.ch/
  //
  TString struttura("");

  struttura += "@md1.pressione(),";
  struttura += "@md1.temperatura(),";
  struttura += "@md1.umidita(),";
  struttura += "@md1.precipitazione()";

  struttura += ",";

  struttura += "@md2.pressione(),";
  struttura += "@md2.temperatura(),";
  struttura += "@md2.umidita(),";
  struttura += "@md2.precipitazione()";

  struttura += ",";

  struttura += "@md3.pressione(),";
  struttura += "@md3.temperatura(),";
  struttura += "@md3.umidita(),";
  struttura += "@md3.precipitazione()";

  struttura += ",";

  struttura += "@md4.pressione(),";
  struttura += "@md4.temperatura(),";
  struttura += "@md4.umidita(),";
  struttura += "@md4.precipitazione()";

  struttura += ":";
  struttura += neuro1;
  if (neuro2>0) { // due strati nascosti
    struttura+=":";
    struttura+=neuro2;
  }
  struttura+=":";

#ifdef PRED_PIOGGIA
  struttura+="@md.precipitazione()";
#else
  struttura+="@md.temperatura()";
#endif

  return struttura;
}

//
// addestra la rete neurale e salva i pesi ottimizzati
// - nAlgo	algoritmo di addestramento, numero da 1 a 6
// - neuro1	numero di neuroni nel primo strato nascosto
// - neuro2	no. di neuroni nel secondo strato nascosto (puo' essere zero)
// - ntrain	no. di cicli di addestramento
// - name	base dei nomi dei file con le n-uple di Root
//
void reteneurale(Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain, TString name="dati")
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
  }
  
  //
  // apertura del file e del suo tree
  // questo sara' il file usato per il training della rete
  //
  TFile rfile(name + ".root");
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)rfile.Get("T");

  //
  // apriamo i file di Root sfasati di 1-4 giorni e usiamoli come file
  // "amici" del precedente, per poter usare nella rete neurale i dati
  // dei giorni precedenti
  //
  tree->AddFriend("T", name + "-1.root");
  tree->AddFriend("T", name + "-2.root");
  tree->AddFriend("T", name + "-3.root");
  tree->AddFriend("T", name + "-4.root");

  //--- creiamo un oggetto della classe "rete neurale" definito da "struttura"
  TString struttura =  nnStruct(neuro1, neuro2);
  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(struttura, tree);

  //
  // collect list of events for the train and test set
  // the neural network is optimized on the train set and
  // comparing the accuracy on the test set one controls for overtraining
  //
  TEventList train;
  TEventList test;
  for (int i = 0; i < tree->GetEntriesFast(); i++) {
    if (i%10 < 5)
      train.Enter(i);
    else
      test.Enter(i);
  }
  train.Print();
  test.Print();

  mlp->SetTrainingDataSet(&train);
  mlp->SetTestDataSet(&test);

  //--- algoritmo di addestramento:
  if (nAlgo == 1) mlp->SetLearningMethod(TMultiLayerPerceptron::kStochastic);
  if (nAlgo == 2) mlp->SetLearningMethod(TMultiLayerPerceptron::kBatch);
  if (nAlgo == 3) mlp->SetLearningMethod(TMultiLayerPerceptron::kSteepestDescent);
  if (nAlgo == 4) mlp->SetLearningMethod(TMultiLayerPerceptron::kRibierePolak);
  if (nAlgo == 5) mlp->SetLearningMethod(TMultiLayerPerceptron::kFletcherReeves);
  if (nAlgo == 6) mlp->SetLearningMethod(TMultiLayerPerceptron::kBFGS);
  if (nAlgo<1 || nAlgo>6) cout << "Algoritmo non valido." << endl;
  cout << "Algoritmo scelto: " << nomeAlgo(nAlgo) << endl;

  //--- addestramento:
  mlp->Train(ntrain, "text, graph, update=10");

  TCanvas* c2 = new TCanvas("c2", "Network analysis");
  c2->SetWindowPosition(50,50);

  //--- mostra i pesi delle sinapsi dopo l'apprendimento
  mlp->Draw();

  //--- stampa informazioni con l'analizzatore della rete neurale
  TMLPAnalyzer mlpa(mlp);
  //--- Initialisation
  mlpa.GatherInformations();
  //--- output to the console
  mlpa.CheckNetwork();
  //--- show average deviation as function of desired output value
  // mlpa.DrawTruthDeviation(0);

  TString nomepesi = nomeFilepesi(name, nAlgo, neuro1, neuro2, ntrain);
  mlp->DumpWeights(nomepesi);
  cout << "File dei pesi salvato col nome " << nomepesi << endl;
}

//
// inizializza la rete neurale con i pesi salvati per studiarla
//
void usarete(Int_t nAlgo, Int_t neuro1, Int_t neuro2, Int_t ntrain, TString name="dati")
{
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    gSystem->Load("libMLP");
  }

  //
  // apertura dei file di Root contenenti le n-uple di dati
  //
  TFile rfile(name + ".root");
  // metto "T" perche' e' il nome del tree contenuto nel file
  TTree *tree = (TTree*)rfile.Get("T");

  //
  // apriamo i file di Root sfasati di 1-4 giorni e usiamoli come file
  // "amici" del precedente, per poter usare nella rete neurale i dati
  // dei giorni precedenti
  //
  tree->AddFriend("T", name + "-1.root");
  tree->AddFriend("T", name + "-2.root");
  tree->AddFriend("T", name + "-3.root");
  tree->AddFriend("T", name + "-4.root");

  meteoData* mdp = 0;
  meteoData* mdp1 = 0;
  meteoData* mdp2 = 0;
  meteoData* mdp3 = 0;
  meteoData* mdp4 = 0;

  tree->SetBranchAddress("md", &mdp);
  tree->SetBranchAddress("md1", &mdp1);
  tree->SetBranchAddress("md2", &mdp2);
  tree->SetBranchAddress("md3", &mdp3);
  tree->SetBranchAddress("md4", &mdp4);

  // definiamo la struttura della rete neurale (dev'essere identica a
  // quella usata in reteneurale())
  TString struttura =  nnStruct(neuro1, neuro2);
  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(struttura, tree);

  //--- carichiamo i pesi delle sinapsi salvati in precedenza
  mlp->LoadWeights(nomeFilepesi(name, nAlgo, neuro1, neuro2, ntrain));

  // loop sul numero di entries nel sample
  Double_t inputs[4*4];
  Double_t output;
  Double_t naive;
  Double_t target;

#ifdef PRED_PIOGGIA
  const double histoMin(-2);
  const double histoMax(2);
  const int histoBins(10);
  const double histoOffs((histoMax-histoMin)/10/2);

  TH1F *hnn = new TH1F("hnn", "T_nn - T", histoBins, histoMin-histoOffs, histoMax-histoOffs);
  TH1F *hnaive = new TH1F("hnaive", "T_naive - T", histoBins, histoMin-histoOffs, histoMax-histoOffs);
#else
  TH1F *hnn = new TH1F("hnn", "T_nn - T", 100, -20, 20);
  TH1F *hnaive = new TH1F("hnaive", "T_naive - T", 100, -20, 20);
#endif

  Int_t pioveNN=0;
  Int_t pioveR=0;
  Int_t falsiPos=0;
  Int_t falsiNeg=0;

  //--- soglia sopra la quale consideriamo piovosa la predizione della rete neurale
  Double_t soglia(0.5);

  for (int i = 0; i < tree->GetEntriesFast(); i++) {
    tree->GetEntry(i);

    if (i%10 < 5) continue;
    //--- test neural network on the test sample

    Int_t ii(0);

    inputs[ii++] = mdp1->pressione();
    inputs[ii++] = mdp1->temperatura();
    inputs[ii++] = mdp1->umidita();
    inputs[ii++] = mdp1->precipitazione();
    
    inputs[ii++] = mdp2->pressione();
    inputs[ii++] = mdp2->temperatura();
    inputs[ii++] = mdp2->umidita();
    inputs[ii++] = mdp2->precipitazione();

    inputs[ii++] = mdp3->pressione();
    inputs[ii++] = mdp3->temperatura();
    inputs[ii++] = mdp3->umidita();
    inputs[ii++] = mdp3->precipitazione();

    inputs[ii++] = mdp4->pressione();
    inputs[ii++] = mdp4->temperatura();
    inputs[ii++] = mdp4->umidita();
    inputs[ii++] = mdp4->precipitazione();

    //--- predizione della rete neurale
    output = mlp->Evaluate(0,inputs);

    //
    // the output is normalized:
    // - the mean of the events to be predicted is subtracted
    // - and the result is divided by the RMS of the same events
    // to compare the net output with the target value, it must be
    // un-normalized using the data that were saved and used for the
    // network initialization
    //
    const Double_t* rmsAndNorm = ((TNeuron*)mlp->fLastLayer[0])->GetNormalisation();
    output = output*rmsAndNorm[0]+rmsAndNorm[1];

#ifdef PRED_PIOGGIA
    //--- la rete prevede pioggia se l'uscita supera 0.5
    output = output > 0.5 ? 1 : 0;
    //--- predizione naive
    naive = mdp1->precipitazione();
    //--- dato metereologico che si vorrebbe prevedere
    target = mdp->precipitazione();
#else
    //--- predizione naive
    naive = mdp1->temperatura();
    //--- dato metereologico che si vorrebbe prevedere
    target = mdp->temperatura();
#endif

    //
    // accumula istogrammi dello scostamento delle predizioni
    //
    hnn->Fill(output - target);
    hnaive->Fill(naive - target);

    if (output > soglia) pioveNN++;
    if (target > soglia) pioveR++;
    if (output > soglia && target <= soglia) falsiPos++;
    if (output <= soglia && target > soglia) falsiNeg++;
  }

  // visualizza confronto tra output e target
  // gStyle->SetOptStat(0000);
  TCanvas* c3 = new TCanvas("c3","Predizione rete neurale e naive", 75,75, 600, 700); 
  c3->Divide(1,2);
  c3->cd(1);
  hnn->DrawCopy();
  c3->cd(2);
  hnaive->DrawCopy();

#if 0
  //
  // vecchio codice salvato
  //

  h1->SetLineWidth(2);
  h2->SetLineWidth(2);
  hnn->SetLineColor(kBlack);
  hnaive->SetLineColor(kRed);

  // legenda
  Double_t ymax = TMath::Max(h1->GetMaximum(),h2->GetMaximum());
  TLegend *legend =
    new TLegend(7.5,0.75*ymax,9.5,0.95*ymax,"Precipitazioni","br");
  legend->AddEntry(h1,"previsione","l");
  legend->AddEntry(h2,"realta`","l");
  legend->Draw();
#endif

#ifdef PRED_PIOGGIA
  cout
    << "Ha piovuto in " << pioveR
    << " giorni contro " << pioveNN
    << " previsti dalla rete."
    << endl
    << "Falsi positivi: " << falsiPos
    << ", falsi negativi: " << falsiNeg
    << endl;
#endif

}

void nnMeteo()
{
  gROOT->Reset(); // fa il reset di tutte le variabili

  // 
  // converte il file testo in un file di ROOT; la variabile numerica
  // e' il numero di giorni descritti dal file se avete piu' di un
  // file di input, ripetete il rigo sopra piu' volte mettendo ogni
  // volta il nome del file (senza ".txt")
  //
  converti("dati", 365);
  
  //--- grafico delle distribuzioni di tutte le variabili
  visualizza("dati");

  //
  // addestra la rete
  // - argomento #1: numero di algoritmo da 1 a 6
  //     1) Stochastic
  //     2) Batch
  //     3) SteepestDescent
  //     4) RibierePolak
  //     5) FletcherReeves
  //     6) BFGS
  // - argomento #2 e #3: neuroni nel 1o e 2o strato nascosto
  // - argomento #4: cicli di apprendimento
  //
  reteneurale(5, 2, 0, 100);

  //
  // ricrea la rete con i pesi dati dall'addestramento, e la mette a
  // confronto con la realta'
  // - argomento #1: numero di algoritmo da 1 a 6
  // - argomento #2 e #3: neuroni nel 1o e 2o strato nascosto
  // - argomento #4: cicli di apprendimento
  //
  usarete(5, 2, 0, 100);
}
